Reverse Monte Carlo modeling of amorphous silicon 
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An implementation of the Reverse Monte Carlo algorithm is presented for the study of amorphous 
tetrahedral semiconductors. By taking into account a number of constraints that describe the 
tetrahedral bonding geometry along with the radial distribution function, we construct a model 
of amorphous silicon using the reverse monte carlo technique. Starting from a completely random 
configuration, we generate a model of amorphous silicon containing 500 atoms closely reproducing 
the experimental static structure factor and bond angle distribution and in improved agreement 
with electronic properties. Comparison is made to existing Reverse Monte Carlo models, and the 
importance of suitable constraints beside experimental data is stressed. 
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I. INTRODUCTION 

The structure of amorphous semiconductors is well rep- 
resented by continuous random network (CRN) model in- 
troduced by Zachariaseni 70 years ago. The CRN model 
has the simplicity that each of the atoms should sat- 
isfy its local bonding requirements and should have as 
small strain as possible in the network which is generally 
characterized by having a narrow bond angle as well as 
bond length distribution. In spite of its apparent sim- 
plicity, the structural modeling of high quality tetrahe- 
dral amorphous semiconductors appears to be quite dif- 
ficult. There have been many models of amorphous sil- 
icon2i2i2i£i§iL2i2ii2iii proposed in the last 30 years which 
include from very simple hand-built model of Polk 2 , com- 
puter generated periodic network model of Guttman 3 
to the complex model of Wooten, Winer and Weaire 
(WWW)^, but most of these models have some limita- 
tions in one way or another in describing the true na- 
ture of the amorphous state. The last method, the so 
called sillium approach of WWW, is based on the strat- 
egy of randomizing and relaxing the network is so far the 
most successful method of producing minimally strained 
CRN. The algorithm, in its modified form developed by 
Djordjevic et. al«2i and Barkema and Mousseau 5 -, can pro- 
duce a CRN which is comparable to experimental results 
and is capable of producing a clean band gap without 
any defect states in the gap. 

In this paper we develop a different approach to model 
amorphous semiconductors known as reverse monte carlo 
(RMC) simulation 12 i 13 i 14 i 15 i 16 i 17 i 18 . Our primary objec- 
tive is to produce structural configurations that are con- 
sistent with experimental data but at the same time we 
go one step further to generate realistic configurations 
for comparison with models obtained via other routes. 
We emphasize that producing realistic models (mean- 
ing models which agree with all experiments) requires 
more than spatial pair correlations, and identify addi- 
tional constraints which lead to realistic models. 

The existing RMC models of amorphous semiconduc- 
tors are found to be inadequate and fail to produce some 
of the basic experimental features of amorphous tetrahe- 



dral semiconductors. Gereben & PusztaiiS*^ have car- 
ried out RMC simulation of tetrahedral semiconductors 
using a number of models ranging from completely dis- 
ordered configuration to randomized diamond structure. 
Although a certain degree of tetrahedral character in the 
bond angle distribution was reflected in their work, most 
of the models show an unphysical peak in bond angle 
distribution around 60° . The work of Walters and New- 
porti* on amorphous germanium made some progress 
toward getting the correct bond angle distribution, but 
the number of 3-fold coordinated atoms are quite high 
in their model and in absence of any discussion on local 
strain and electronic properties it is difficult to say how 
reliable their models are when it comes looking at the 
electronic properties. 



A developing area where RMC may be applied suc- 
cessfully is for modeling amorphous materials exhibiting 
medium-range order (MRO). Such MRO is characterized 
by the existence of 10 — 20A scale structure. Recent de- 
velopments in fluctuation electron microscopy (FEM)iS 
and its application on amorphous germanium and silicon 
have indicated that computer generated CRN model of 
these materials lack the characteristic signature of MRO. 
Since RMC is based on experimental data, it provides 
a promising scheme to model amorphous materials hav- 
ing medium-range order by including the experimentally 
measured FEM signal as input data to augment pair cor- 
relations. 



The plan of the paper is as follows. In section [HI 
we briefly mention the basic philosophy of reverse monte 
carlo modeling and some of its salient features. This is 
followed by role of constraints in RMC modeling in sec- 
tion lHII where we illustrate how a set of judiciously chosen 
constraints can be used to construct a reliable model of 
amorphous silicon. Finally we compare our results with 
those obtained from earlier RMC models and a model 
obtained via WWW algorithm. 
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II. BASICS OF RMC 

The RMC method has been described in detail else- 
where^. Here we briefly outline the basic philosophy 
of RMC. At the very basic level, RMC is a technique 
for generating structural configurations based on experi- 
mental data. The logic is very appealing: any model of 
a complex material worthy of further study should, at 
a minimum, agree with what is known (that is, the ex- 
periments). By construction, the RMC scheme enforces 
this (and for contrast, a molecular dynamics simulation 
may not). In an ideal implementation, one should find 
a model agreeing with all known information, but this is 
not easy to accomplish, though we make some progress 
below. The approach was originally developed by Mc- 
Greevy & Puszta i 12 i 15 for liquid and glassy materials for 
lack of different routes to explore experimental data but 
in recent years progress has been made toward modeling 
crystalline systems as well2£. Starting with a suitable 
configuration, atoms are displaced randomly using the 
periodic boundary condition until the input experimental 
data (either the structure factor or the radial distribution 
function) match with the data obtained from the gener- 
ated configuration. This is achieved by minimizing a cost 
function which consists of either structure factor or radial 
distribution function along with some appropriately cho- 
sen constraints to restrict the search space. Consider a 
system having N number of atoms with periodic bound- 
ary condition. One can construct a generalized cost func- 
tion for an arbitrary configuration by writing : 

K M L 
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where r\\ is related to the uncertainty associated with 
the determination of experimental data points as well as 
the relative weight factor for each set of different exper- 
imental data. The quantity Q is the appropriate gener- 
alized variable associated with experimental data F(Q) 
and Pi is the penalty function associated with each con- 
straint. For example, in case of radial distribution func- 
tion and structure factor, Q has the dimension of length 
and inverse length respectively. In order to avoid the 
atoms getting too close to each other, a certain cut-off 
distance is also imposed which is typically of the order 
of interatomic spacing. In RMC modeling, this is usu- 
ally obtained from the radial distribution function by 
Fourier transform of the measured structure factor. This 
is equivalent to adding a hard sphere potential cut-off in 
the system which prevents the catastrophic build up of 
potential energy. 

In spite of the fact that RMC has been applied to many 
different types of systems - liquid, glasses, polymer and 
magnetic materials, questions are often raised about the 
reliability of results obtained from RMC simulation. The 
method has never been accepted without some degree of 
controversy and the most popular criticism is the lack of 
unique solution from RMC. RMC can produce multiple 



configurations having the same pair correlation function. 
This lack of uniqueness, however, is not surprising, since 
usually only the pair correlation function or structure 
factor is used in modeling the structure, while there ex- 
ists an infinite hierarchy of higher order correlation func- 
tions carrying independent structural information are ne- 
glected. In other words, RMC samples from the space of 
all models consistent with some limited body of data - in 
its simplest form (analyzing a single experiment) RMC 
is an ideal gauge of how non-specific the data is with re- 
spect to identification of an atomistic model. If the mod- 
eler possesses a priori information independent of that 
implicit in the experiment being fit to, it is necessary to 
add this information to the modeling in some fashion to 
receive a model in joint agreement with the experiment 
and the additional information. 



III. A NEW RMC MODEL 

We begin by including the minimal information that 
is necessary to model a configuration of a-Si. In so do- 
ing, we use the radial distribution function obtained from 
a high quality model of amorphous silicon. This latter 
model was generated by Barkema and Mousseau* using 
a modified form of WWW algorithm 4 having bond angle 
distribution close to 10° with 100% 4-fold coordination. 
In addition to this RDF, we also impose the conditions 
that the average bond angle of all the triplets Si-Si-Si 
should be near 109.5° and the corresponding root mean 
square deviation should be no less than 10°. The num- 
ber of 4-fold coordinated atoms is driven to a specified 
value during the simulation by including a constraint on 
the average coordination number. It is to be noted that 
while there is no limit to the number of constraints that 
can be included in the system, there is no guarantee that 
mere inclusion of more constraints will necessarily give 
better results. Forcing a completely random configura- 
tion with too many competing constraints may cause the 
configuration to be trapped in the local minimum of the 
function £ and may prevent the system from exploring 
a large part of the search space. By adding only the 
essential constraints that describe the chemical and ge- 
ometrical nature of the bonding correctly, Eq.^ can be 
written as : 



M 

i=i 

+ A 2 (0 o -fl) 2 

+ \ 3 (S9 -S9) 2 

+ A 4 {1 - Q(x - x c )} 

+ A 5 {0 o -0} 2 (2) 

where, 
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FIG. 1: Structure factor obtained from a RMC (+) model 
containing 500 atoms of a-Si. The solid line is obtained from a 
WWW sample of identical size and number density of atoms. 
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FIG. 2: Structure factor obtained for a 500-atom model of 
a-Si from RMC (solid line) and the experiment of Laaziri et 
al. as indicated in the figure. 



56 = V((0-0 ) 2 ) 

In Eq.0 and 89 are the average angle and the rms 
deviation while <fi an d </>o are the current and proposed 
concentration of the 4- fold coordinated atoms. It is im- 
portant to note that each of the terms in Eq.|21 is non- 
negative, and should decrease ideally to zero during the 
course of minimization. Since the cost of energy associ- 
ated with the bond length relaxation is more than the 
bond angle relaxation, atomic arrangements with large 
bond angle distribution but having correct RDF fre- 
quently result. The coefficients, Ai to A3, for the different 
terms in Eq.|21 can be chosen appropriately to minimize 
this effect. In general the coefficients A are constant dur- 
ing the course of simulation but the minimization pro- 
cedure can be slightly accelerated by making them vary 
in such a way that the contribution from each of the 
term are of the same order during the course of sim- 
ulation. The coefficient A4 is usually assigned a large 
value in order to include a hard sphere cut-off as men- 
tioned earlier so that no two particles can come closer 
to x c while the coefficient A5 maintains the number of 
4-fold coordinated atoms to a specified value. In RMC 
simulation of amorphous tetrahedral semiconductors one 
usually encounters the problem of having a pronounced 
peak at 60° . This peak is a characteristic feature of un- 
constrained RMC simulation and is due to the formation 



of equilateral triangles by three atoms. In the work of 
Gereben and PusztaiiSiii, attempts were made to over- 
come this difficulty by constraining the bond angle dis- 
tribution as well as by making an initial configuration 
which is 100% 4-fold obtained from a diamond lattice. 
The resulting structure is, however, found to be unstable 
and on relaxation using a suitable potential, the config- 
uration tends to get back toward the starting structure, 
i.e., randomized diamond in this case—. In the approach 
of Walter and Newport 14 , the initial random configura- 
tion was examined and any "triples" , i.e. , three atoms 
forming an equilateral triangle was removed before the 
beginning of RMC fit. By selective removal of such un- 
wanted triplets, they have been able to generate con- 
figuration of a-Ge without having a peak at 60°. The 
approach that we have taken in our work is more general 
and starts with a completely random configuration. This 
eliminates, in the first place, any possible local ordering 
that may exist in the starting structure (e.g., randomized 
diamond structure retains the memory of tetrahedral or- 
dering). Furthermore, we have not included or excluded 
any special configuration in our starting structure, e.g., 
three atoms forming an equilateral triangle. Based on ex- 
perimental consideration, we have included only the key 
features of amorphous tetrahedral semiconductors - an 
average bond angle of 109.5° having rms deviation of 10° 
which is consistent with the RDF obtained from a WWW 
relaxed model used in our calculation. For the 500-atom 
model reported in this work, we have chosen a cubic box 
of length 21.18A which corresponds to number density 
0.0526 atom/A 3 . The initial configuration is generated 
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FIG. 3: The bond angle distribution functions (BADF) 
for 500-atom a-Si from constrained RMC (dashed line) and 
WWW model (solid line). The rms deviation for the models 
are 12.5° and 9.9° respectively. 

randomly so that no two atoms can come closer to 2.0A. 
The configuration is then relaxed by moving the atoms 
to minimize the cost function £. In addition to applying 
standard monte carlo moves in which a single or a group 
of atoms is randomly displaced, a variety of monte carlo 
moves have been implemented in our work. For example, 
in one of such moves, a 3-fold or 5-fold atom is selected 
and the nearest neighbor distance is examined. If the 
distance is greater than 2.7A, the neighboring atom is 
displaced in order to bring the distance within a radius 
of 2.7k. The maximum displacement of a monte carlo 
move is limited to 0.2-0.4A throughout the simulation. 
Since we are interested in the electronic structure as well, 
we confine ourselves within a reasonable system size for 
studying the generated structure using a first principles 
density functional Hamiltonian. The density functional 
calculations were performed within the local density ap- 
proximation (LDA) using the local basis first principles 
code Siesta^. We have used a non self-consistent ver- 
sion of density functional theory based on the lineariza- 
tion of the Kohn-Sham equation by Harris functional ap- 
proximation 23 along with the parameterization of Perdew 
and Zunger— for the exchange-correlation functional. 



IV. RESULTS 

The results for the model including all the constraints 
are presented in Figs. 11141 Since the structure factor is 
generally considered to be more sensitive to an arbitrary 
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FIG. 4: The electronic density of states (EDOS) of 500-atom 
model of a-Si obtained from RMC simulation described in the 
text. The Fermi level is at E=0. 



small change in the atomic positions than the radial dis- 
tribution function, we have plotted the structure factors 
for the constrained RMC and WWW model in Fig-QJ It 
is evident from Fig.^ that the agreement between the 
RMC and WWW model is very good both for small and 
large values of Q. In order to further justify the credibil- 
ity of our model, we have plotted in FigJSJthe structure 
factor from the experiment of Laaziri et. al 25 along with 
the same obtained from our RMC model. Once again we 
find that the agreement between the structure factor from 
RMC and the experimental results is quite good except 
for the few points near the first peak. It is very tempting 
to think this deviation as a finite size effect coming from 
the finiteness of our model. We have therefore calculated 
the structure factor for WWW models containing 300 to 
4096 atoms of Si but the deviation continues to remain. 
Holender and Morganii also observed similar deviation 
near the first peak in their work with a much larger model 
containing 13824 atoms which was compared with the ex- 
perimental data obtained by Fortner and Lannin2&. 

In Fig.|3J we have plotted the bond angle distributions 
(BADF) for both the RMC and WWW model. As we 
have discussed in section II, the radial distribution func- 
tion or structure factor can not alone provide all the 
necessary information that are needed to characterize an 
atomic configuration obtained from a reverse monte carlo 
simulation. A further characterization beyond pair cor- 
relation function is therefore vital and necessitates the 
need for getting some idea about the 3-body correlation 
function. It is clear from the Fig .|3 that the distribution 
obtained from the RMC model follows the tetrahedral 
character observed in amorphous semiconductors. The 
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average bond angle in this case is found to be 109.01° 
with rms deviation of 12.5°. An important aspect of the 
bond angle distribution in Figgis that most of the an- 
gles are lying between 70°-150° compared to 80°-140° in 
WWW case. We emphasize at this point that the ear- 
lier works on modeling amorphous tetrahedral semicon- 
ductors using RMC predicted a much wider bond angle 
distributions. Gereben and Pusztaii^ have observed a 
pronounced, unphysical peak at 60° except for the model 
starting with diamond structure while Walters and New- 
por1jA4 have reported a bond angle distribution of a-Gc 
which is as wide as 60°-180°. It is an important develop- 
ment here that by adding three more constraints (A2, A3, 
and A4) we have achieved a significantly improved results. 
Both the radial and the bond angle distribution func- 
tions reported here are at par with the results obtained 
from molecular dynamics simulation and is comparable 
to those obtained from WWW model. The fact that the 
inclusion of these two constraints leads to a significant 
improvement is not surprising. For a large continuous 
random network (CRN) model of amorphous tetrahe- 
dral semiconductor, one can approximate the bond angle 
distribution as nearly GaussianSl. This approximated 
Gaussian distribution can defined by the first two mo- 
ments of the distribution function. By specifying these 
two moments as constraints in Eq.|21 we correctly de- 
scribe the tetrahedral bonding geometry of the atoms 
which along with the radial distribution function pro- 
duces a configuration more realistic than those obtained 
from models based on RDF or structure factor only. This 
suggests that in addition to the radial distribution of the 
atoms, one needs to include some relevant information 
about the nature of 3-body correlation among the atoms 
to construct a realistic configuration. 

Having studied the radial and bond angle distribution 
we now address the electronic density of states calcula- 
tions. While the width of the bond angle distribution 
function (BADF) and the structure factor together in- 
deed gives some idea about the quality of the model, 
some of the features e.g., the existence of spectral gaps 
and the position of defects states in the spectrum can be 
studied by looking at the electronic density of states only. 
The structure obtained from RMC simulation is first re- 
laxed using the density functional code Siesta and is 
found to be close to an energy minimum in the local den- 
sity approximation (LDA). This is an important test for 
determining the stability of the structure obtained from 
RMC simulation and as far as we are concerned almost 
all earlier works on RMC have completely neglected this 
issue. In Fig. 01 we have plotted the electronic density 
of states (EDOS) for the constrained model. The EDOS 
appears with all the characteristic features of a-Si with 
the exception of a clean gap in the spectrum. This be- 
havior is not unexpected in view of the fact that 88% 
of the total atoms are found to be 4-fold coordinated 
with an average coordination number 3.85. The pres- 
ence of the defect states makes the gap noisy and at the 



same time the use of LDA underestimates the size of 
the gap. This EDOS is in significantly better agreement 
with optical measurements than conventional RMC mod- 
els with much higher defect concentrations and spurious 
bond angles. It is interesting to observe that the average 
coordination number from our model is very close to the 
experimental value of 3.88 reported by Laaziri et. al~. 



V. CONCLUSIONS 

We have presented a model of amorphous silicon based 
on reverse monte carlo simulation. One of the novel fea- 
tures of our model is to start with a completely random 
structure and then to relax toward a realistic configu- 
ration by adding a number of physically relevant con- 
straints. The characteristic features of the tetrahedral 
bonding are taken into account by adding constraints on 
average bond angle and its deviation from the mean while 
the number of 4-fold coordinated atoms is maintained at 
a specified value by further use of a constraint on aver- 
age coordination number. The radial and the bond angle 
distribution obtained from our model is found to be in 
excellent agreement with a high quality CRN model pro- 
duced by WWW algorithm. We have also compared the 
structure factor with the experimental data obtained by 
Laaziri et al. and observed a reasonably good agreement. 
By relaxing the model using the first principles density 
function code Siesta, we find that the model is close 
to the energy minimum for LDA and is stable. The elec- 
tronic density of states (EDOS) obtained from our model 
contains all the essential feature of amorphous silicon in- 
cluding a signature of the band gap. Although the model 
does not produce a clean gap in the spectrum, the quality 
of the EDOS is at par with models obtained from molec- 
ular dynamics simulation. Our RMC algorithm presents 
a significant improvement on previous RMC studies and 
makes it possible to compare for the first time, albeit 
qualitatively, the structural and electronic properties of 
RMC models with its WWW counterpart. We expect 
that further developments toward this direction will even- 
tually make RMC as an useful modeling tool incorporat- 
ing experimental information and can be used effectively 
without any criticisms in modeling complex materials. 
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